pacman::p_load(sf, tmap, tidyverse, rvest)Thematic Mapping and GeoVisualisation with R
Overview
Plotting functional and truthful choropleth maps by using an R package called tmap package.
Getting Started
Install and launching R packages
The code chunk below uses p_load() of pacman package to check if the following packages are installed in the computer. If they are, then they will be launched into R.
sf - Provides the core tools for handling spatial data
tmap - A package used for producing maps for visualisation
readr - A package for importing delimited text files
tidyr - A package for tidying and reshaping data into a clean format
dplyr - A package for wrangling and manipulating data
rvest - A web scraping package to download and parse data from websites
To Note: readr, tidyr and dplyr are part of tidyverse package so there is no need to load these packages individually.
Importing the Data
Datasets
The following datasets will be used:
Master Plan 2019 Subzone Boundary (No Sea) kml data file (Master plan 2019)
Singapore Residents by Planning Area / Subzone, Age Group, Sex and Type of Dwelling, June 2024 csv file (respopagesextod2024)
Importing Geospatial Data into R
st_read() function of the sf package will be used to import MP14_SUBZONE_WEB_PL shapefile into mpsz.
mpsz <- st_read("data/geospatial/MasterPlan2019SubzoneBoundaryNoSeaKML.kml")Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source
`C:\Users\zongy\OneDrive\Desktop\SMU\ISSS626 - Geospatial Analytics\zongyin-tan\ISSS626-Geospatial-zytan\Hands-on_Ex\Hands-on_Ex01b\data\geospatial\MasterPlan2019SubzoneBoundaryNoSeaKML.kml'
using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84Tidying the Data
A function called extract_kml_field is created to extract values such as REGION_N, PLN_AREA_N, SUBZONE_N, and SUBZONE_C from the HTML Description field. This is done using the code chunk below.
extract_kml_field <- function(html_text, field_name) {
if (is.na(html_text) || html_text == "") return(NA_character_)
page <- read_html(html_text)
rows <- page %>% html_elements("tr")
value <- rows %>%
keep(~ html_text2(html_element(.x, "th")) == field_name) %>%
html_element("td") %>%
html_text2()
if (length(value) == 0) NA_character_ else value
}The code chunk below then applies this function to create new columns REGION_N, PLN_AREA_N, SUBZONE_N, and SUBZONE_C in the dataset. The raw Name and Description fields are removed, and geometry is moved to the last column for better structure.
mpsz <- mpsz %>%
mutate(
REGION_N = map_chr(Description, extract_kml_field, "REGION_N"),
PLN_AREA_N = map_chr(Description, extract_kml_field, "PLN_AREA_N"),
SUBZONE_N = map_chr(Description, extract_kml_field, "SUBZONE_N"),
SUBZONE_C = map_chr(Description, extract_kml_field, "SUBZONE_C")
) %>%
select(-Name, -Description) %>%
relocate(geometry, .after = last_col())We can view the content using the following code chunk.
mpszSimple feature collection with 332 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84
First 10 features:
REGION_N PLN_AREA_N SUBZONE_N SUBZONE_C
1 CENTRAL REGION BUKIT MERAH DEPOT ROAD BMSZ12
2 CENTRAL REGION BUKIT MERAH BUKIT MERAH BMSZ02
3 CENTRAL REGION OUTRAM CHINATOWN OTSZ03
4 CENTRAL REGION DOWNTOWN CORE PHILLIP DTSZ04
5 CENTRAL REGION DOWNTOWN CORE RAFFLES PLACE DTSZ05
6 CENTRAL REGION OUTRAM CHINA SQUARE OTSZ04
7 CENTRAL REGION BUKIT MERAH TIONG BAHRU BMSZ10
8 CENTRAL REGION DOWNTOWN CORE BAYFRONT SUBZONE DTSZ12
9 CENTRAL REGION BUKIT MERAH TIONG BAHRU STATION BMSZ04
10 CENTRAL REGION DOWNTOWN CORE CLIFFORD PIER DTSZ06
geometry
1 MULTIPOLYGON (((103.8145 1....
2 MULTIPOLYGON (((103.8221 1....
3 MULTIPOLYGON (((103.8438 1....
4 MULTIPOLYGON (((103.8496 1....
5 MULTIPOLYGON (((103.8525 1....
6 MULTIPOLYGON (((103.8486 1....
7 MULTIPOLYGON (((103.8311 1....
8 MULTIPOLYGON (((103.8589 1....
9 MULTIPOLYGON (((103.8283 1....
10 MULTIPOLYGON (((103.8552 1....Importing Attribute Data into R
we will import respopagesextod2024.csv file into RStudio and save the file into an tibble dataframe called popdata.
The task will be performed by using read_csv() function of readr package using the following code chunk.
popdata <- read_csv("data/aspatial/respopagesextod2024.csv")Rows: 100928 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): PA, SZ, AG, Sex, TOD
dbl (2): Pop, Time
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Data Preparation
Before creating the thematic map, we need to prepare a data table with year 2020 values. The data table should include following variables:
YOUNG - age group 0 to 4 until age group 20 to 24,
ECONOMY ACTIVE - age group 25-29 until age group 60-64,
AGED - age group 65 and above,
TOTAL - all age group, and
DEPENDENCY - the ratio between young and aged against economy active group
Data Wrangling
The following data wrangling and transformation functions will be used:
pivot_wider() of tidyr package, and
mutate(), filter(), group_by() and select() of dplyr package
popdata2024 <- popdata %>%
group_by(PA, SZ, AG) %>%
summarise(`POP` = sum(`Pop`)) %>%
ungroup()%>%
pivot_wider(names_from=AG,
values_from=POP) %>%
mutate(YOUNG = rowSums(.[3:6])
+rowSums(.[12])) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[7:11])+
rowSums(.[13:15]))%>%
mutate(`AGED`=rowSums(.[16:21])) %>%
mutate(`TOTAL`=rowSums(.[3:21])) %>%
mutate(`DEPENDENCY` = (`YOUNG` + `AGED`)
/`ECONOMY ACTIVE`) %>%
select(`PA`, `SZ`, `YOUNG`,
`ECONOMY ACTIVE`, `AGED`,
`TOTAL`, `DEPENDENCY`)`summarise()` has grouped output by 'PA', 'SZ'. You can override using the
`.groups` argument.Joining the attribute data and geospatial data
Before performing the georelational join, we need to standardise the text format of the PA and SZ fields by converting all values to uppercase. This step is necessary because the current values contain a mix of upper- and lowercase letters, whereas the SUBZONE_N and PLN_AREA_N fields are already stored entirely in uppercase.
popdata2024 <- popdata2024 %>%
mutate_at(.vars = vars(PA, SZ),
.funs = list(toupper)) %>%
filter(`ECONOMY ACTIVE` > 0)left_join() of dplyr is used to join the geographical data and attribute table using planning subzone name e.g. SUBZONE_N and SZ as the common identifier.
mpsz_pop2024 <- left_join(mpsz, popdata2024,
by = c("SUBZONE_N" = "SZ"))Finally, writing after data preparation is completed into an rds file.
write_rds(mpsz_pop2024, "data/rds/mpsz_pop2024.rds")Choropleth Mapping Geospatial Data Using tmap
Choropleth mapping is a technique used to represent enumeration units such as countries, provinces, states or census areas by filling them with patterns or graduated colors.
Two approaches can be used to prepare thematic map using tmap, they are:
Plotting a thematic map quickly by using qtm().
Plotting highly customisable thematic map by using tmap elements.
Plotting a choropleth map quickly by using qtm()
The easiest and quickest to draw a choropleth map using tmap is using qtm(). It is concise and provides a good default visualisation. The following chunk code will be used for that.
tmap_mode("plot")ℹ tmap mode set to "plot".qtm(shp = mpsz_pop2024,
fill = "DEPENDENCY")
From the code above, we learn two key points:
tmap_mode(“plot”) is used to create a static map. If we want an interactive map, we would use tmap_mode(“view”) instead.
The fill argument tells the map which attribute to display. In this case, the DEPENDENCY field.
Creating a choropleth map by using tmap’s elements
While qtm() is useful for creating choropleth maps quickly and easily, its main limitation is the lack of flexibility in controlling the aesthetics of individual layers. To produce a high-quality, publication-ready choropleth map as shown below, it is better to use tmap’s drawing elements, which offer much greater customisation.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
n = 5,
values = "brewer.blues"),
fill.legend = tm_legend(
title = "Dependency ratio")) +
tm_title("Distribution of Dependency Ratio by planning subzone") +
tm_layout(frame = TRUE) +
tm_borders(fill_alpha = 0.5) +
tm_compass(type="8star", size = 2) +
tm_scalebar() +
tm_grid(alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from Urban Redevelopment Authorithy (URA)\n and Population data from Department of Statistics DOS",
position = c("left", "bottom"))
Drawing a base map
In the code chunk below, tm_shape() is used to define the input data (i.e mpsz_pop2024) and tm_polygons() is used to draw the planning subzone polygons.
tm_shape(mpsz_pop2024) +
tm_polygons()
Drawing a choropleth map using tm_polygons()
To draw a choropleth map showing the geographical distribution of a selected variable by planning subzone, we just need to assign the target variable such as Dependency to tm_polygons().
tm_shape(mpsz_pop2024)+
tm_polygons(fill = "DEPENDENCY")
To Note:
The default interval binning used to draw the choropleth map is called “pretty”.
The default colour scheme used is
blues3of ColorBrewer.By default, Missing value will be shaded in grey.
Drawing a choropleth map using tm_fill() and *tm_border()
tm_polygons() is a wraper of tm_fill() and tm_border(). tm_fill() shades the polygons by using the default colour scheme and tm_borders() adds the borders of the polygon features onto the choropleth map.
The code chunk below draws a choropleth map by using tm_fill() alone.
tm_shape(mpsz_pop2024)+
tm_fill("DEPENDENCY")
The planning subzones are shared according to the respective dependecy values but the boundaries are missing. To add the boundary of the planning subzones, tm_borders() will be used as shown in the code chunk below.
tm_shape(mpsz_pop2024)+
tm_fill("DEPENDENCY") +
tm_borders()
Light-gray border lines have been added to the choropleth map using tm_borders(). The fill_alpha argument controls transparency, with values ranging from 0 (fully transparent) to 1 (fully opaque, default). In addition, tm_borders() also allows customisation of border appearance through three arguments: col to set the border colour, lwd to adjust the line width (default is 1), and lty to define the line type (default is “solid”).
tm_shape(mpsz_pop2024)+
tm_fill("DEPENDENCY") +
tm_borders(col = "grey60",
lwd = 0.1,
lty = "dashed")
Data classification methods of tmap
Most choropleth maps employ some methods of data classification. The point of classification is to take a large number of observations and group them into data ranges or classes.
tmap provides a total ten data classification methods, namely: fixed, sd, equal, pretty (default), quantile, kmeans, hclust, bclust, fisher, and jenks.
To define a data classification method, the style argument of tm_fill() or tm_polygons() will be used.
Plotting choropleth maps with built-in classification methods
The code chunk below shows a quantile data classification that used 5 classes.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
n = 5)) +
tm_borders(fill_alpha = 0.5)
In the code chunk below, equal data classification method is used.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "equal",
n = 5)) +
tm_borders(fill_alpha = 0.5)
Notice that the distribution of quantile data classification method are more evenly distributed then equal data classification method.
Plotting choropleth map with custom break
In tmap, category breaks are usually set automatically. If we want to control them, we can use the breaks argument in tm_scale_intervals(). Breaks must include both a minimum and maximum value, so to create n categories, we need to provide n+1 break points in increasing order.
Before deciding on the break points, it’s good practice to check the descriptive statistics of the variable. For example, the code below shows a summary of the DEPENDENCY field:
summary(mpsz_pop2024$DEPENDENCY) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.1905 0.7450 0.8377 0.8738 0.9366 12.7500 94 With reference to the results above, we set break point at 0.60, 0.70, 0.80, and 0.90. In addition, we also need to include a minimum and maximum, which we set at 0 and 100. Our breaks vector is thus (0, 0.60, 0.70, 0.80, 0.90, 1.00)
Now, we will plot the choropleth map by using the code chunk below.
tm_shape(mpsz_pop2024)+
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
breaks = c(0, 0.60, 0.70, 0.80, 0.90, 1.00))) +
tm_borders(fill_alpha = 0.5)Warning: Values have found that are higher than the highest break. They are
assigned to the highest interval
Colour Scheme
tmap supports colour ramps either defined by the user or a set of predefined colour ramps from the RColorBrewer package.
Using ColourBrewer palette
To change the colour, we assign the preferred colour to palette argument of values as shown in the code chunk below.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
n = 5,
values = "brewer.greens")) +
tm_borders(fill_alpha = 0.5)
To reverse the colour shading, we just add a “-” prefix.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
n = 5,
values = "-brewer.greens")) +
tm_borders(fill_alpha = 0.5)
Cartographic Furniture
Beside map style, tmap also also provides arguments to draw other map furniture such as compass, scale bar and grid lines.
In the code chunk below, tm_compass(), tm_scale_bar(), tm_grid() and tm_credit() are used to add compass, scale bar, grid lines and data sources onto the choropleth map.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
n = 5)) +
tm_borders(fill_alpha = 0.5) +
tm_compass(type="8star", size = 2) +
tm_scalebar() +
tm_grid(lwd = 0.1, alpha = 0.2) +
tm_credits("Source: data.gov.sg & singstat",
position = c("left", "bottom"))
Map Layout
A map layout brings together all the elements of a map such as the background, frame, typography, scale, aspect ratio, etc. into a clear and cohesive presentation.
We can refine and customize the layout using the tm_layout() function. In the next 2 sections, we will explore the most commonly used arguments of this function, using the dependency choropleth map as an example.
Map Legend
In tmap, several legend options are provided to change the placement, format and appearance of the legend.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
n = 5),
fill.legend = tm_legend(
title = "Dependency ratio")) +
tm_pos_auto_in() +
tm_borders(fill_alpha = 0.5) +
tm_compass(type="8star", size = 2) +
tm_scalebar() +
tm_grid(lwd = 0.1, alpha = 0.2) +
tm_credits("Source: data.gov.sg & singstat",
position = c("left", "bottom"))
Map style
tmap allows a wide variety of layout settings to be changed. They can be called by using tmap_style().
The code chunk below shows the classic style is used.
tm_shape(mpsz_pop2024) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
n = 5,
values = "-brewer.greens")) +
tm_borders(fill_alpha = 0.5) +
tmap_style("natural")style set to "natural"other available styles are: "white" (tmap default), "gray", "cobalt", "albatross", "beaver", "bw", "classic", "watercolor"tmap v3 styles: "v3" (tmap v3 default), "gray_v3", "natural_v3", "cobalt_v3", "albatross_v3", "beaver_v3", "bw_v3", "classic_v3", "watercolor_v3"
Drawing Small Multiple Choropleth Maps
Small multiple maps, also known as facet maps, consist of several maps arranged side by side or stacked vertically. They are particularly useful for visualizing how spatial relationships change with respect to another variable, such as time.
In tmap, small multiple maps can be plotted in three ways:
by assigning multiple values to at least one of the asthetic arguments,
by creating multiple stand-alone maps with tmap_arrange(), and
by defining a group-by variable in tm_facets().
By assigning multiple values to at least one of the aesthetic arguments
Small multiple choropleth maps are created by assigning two variables to the visual variable using the code chunk below
tm_shape(mpsz_pop2024) +
tm_polygons(
fill = c("YOUNG", "AGED"),
fill.legend =
tm_legend(position = tm_pos_in(
"right", "bottom")),
fill.scale = tm_scale_intervals(
style = "equal",
n = 5,
values = "brewer.blues")) +
tm_borders(fill_alpha = 0.5) +
tmap_style("natural")style set to "natural"other available styles are: "white" (tmap default), "gray", "cobalt", "albatross", "beaver", "bw", "classic", "watercolor"tmap v3 styles: "v3" (tmap v3 default), "gray_v3", "natural_v3", "cobalt_v3", "albatross_v3", "beaver_v3", "bw_v3", "classic_v3", "watercolor_v3"
By arrange multiples choropleth maps in a grid layout
multiple choropleth maps are created and tmap_arrange() is used to arrange them in a grid layout.
youngmap <- tm_shape(mpsz_pop2024)+
tm_polygons(fill = "YOUNG",
fill.legend = tm_legend(
position = tm_pos_in(
"right", "bottom"),
item.height = 0.8),
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.blues")) +
tm_borders(fill_alpha = 0.5) +
tm_title("Distribution of young population")
agedmap <- tm_shape(mpsz_pop2024)+
tm_polygons(fill = "AGED",
fill.legend = tm_legend(
position = tm_pos_in(
"right", "bottom"),
item.height = 0.8),
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.blues")) +
tm_borders(fill_alpha = 0.5) +
tm_title("Distribution of aged population")
tmap_arrange(youngmap, agedmap, asp=1, ncol=2)
By defining a group-by variable in tm_facets()
Multiple small choropleth maps are created by using tm_facets().
tm_shape(mpsz_pop2024) +
tm_fill(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.blues")) +
tm_facets(by = "REGION_N",
nrow = 2,
ncols = 3,
free.coords=TRUE,
drop.units=TRUE) +
tm_layout(legend.show = TRUE,
title.position = c("center", "center"),
title.size = 20) +
tm_borders(fill_alpha = 0.5)
Mapping Spatial Object Meeting a Selection Criterion
Instead of creating small multiple choropleth map, you can also use filter() of dplyr package to select geographical area of interest and plot a choropleth map focus only on the selected region.
mpsz_pop2024 %>%
filter(REGION_N == "CENTRAL REGION") %>%
tm_shape() +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.greens"),
fill.legend = tm_legend()) +
tm_borders(fill_alpha = 0.5)
Complementing Thematic Map with Statistical Chart
Maps and charts work well together because they highlight different aspects of the same data. Maps are good for showing spatial patterns and relationships, while charts make it easy to see numbers, trends, and comparisons. Using both gives a clearer and more engaging view of the data.
In tmap, we can combine maps with statistical charts by using the fill.chart argument and the legend chart feature, as shown in the code chunk below.
mpsz_pop2024 %>%
filter(REGION_N == "CENTRAL REGION") %>%
tm_shape() +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.greens"),
fill.legend = tm_legend(),
fill.chart = tm_chart_box()) +
tm_borders() +
tm_layout(asp = 0.8)
In the code chunk below, We improve the visual representation further by highlighting and lebaling the outliers on the choropleth map.
mpsz_selected <- mpsz_pop2024 %>%
filter(REGION_N == "CENTRAL REGION")
stats <- boxplot.stats(mpsz_selected$DEPENDENCY)
outlier_vals <- stats$out
outlier_sf <- mpsz_selected[mpsz_selected$DEPENDENCY %in% outlier_vals, ]
tm_shape(mpsz_selected) +
tm_polygons(fill = "DEPENDENCY",
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.blues"),
fill.legend = tm_legend(),
fill.chart = tm_chart_box()) +
tm_borders(fill_alpha = 0.5) +
tm_shape(outlier_sf) +
tm_borders(col = "red", lwd = 2) +
tm_text("SUBZONE_N", col = "red", size = 0.7) +
tm_layout(asp = 0.8)
Creating Interactive Map
Interactive maps allow users to explore data by zooming, panning, clicking on locations, and adding overlays, making the experience more dynamic than static maps. With tmap, you can easily switch between static and interactive views using tmap_mode(), depending on your analysis needs.
The code chunks below show how to build an interactive map.
region_selected <- mpsz_pop2024 %>%
filter(REGION_N == "CENTRAL REGION")
region_bbox <- st_bbox(region_selected)
stats <- boxplot.stats(region_selected$DEPENDENCY)
outlier_vals <- stats$out
outlier_sf <- region_selected[region_selected$DEPENDENCY %in% outlier_vals, ]
tmap_mode("view")ℹ tmap mode set to "view".tm_shape(region_selected,
bbox = region_bbox) +
tm_fill("DEPENDENCY",
id = "SUBZONE_N",
popup.vars = c(
"Name" = "SUBZONE_N",
"Dependency" = "DEPENDENCY")) +
tm_borders() +
tm_shape(outlier_sf) +
tm_borders(col = "red", lwd = 2)
<<<<<<< HEAD
=======
>>>>>>> 81336bfc1c4053bed01e1b7bc57e3309e34bd5e4
tmap_mode("plot")ℹ tmap mode set to "plot".The interactive map can be confusing if users zoom in and out too freely. To prevent this, the set_zoom_limits argument is used to restrict how far users can zoom in or out of the map.
region_selected <- mpsz_pop2024 %>%
filter(REGION_N == "CENTRAL REGION")
region_bbox <- st_bbox(region_selected)
stats <- boxplot.stats(region_selected$DEPENDENCY)
outlier_vals <- stats$out
outlier_sf <- region_selected[region_selected$DEPENDENCY %in% outlier_vals, ]
tmap_mode("view")ℹ tmap mode set to "view".tm_shape(region_selected,
bbox = region_bbox) +
tm_fill("DEPENDENCY",
id = "SUBZONE_N",
popup.vars = c(
"Name" = "SUBZONE_N",
"Dependency" = "DEPENDENCY")) +
tm_borders() +
tm_shape(outlier_sf) +
tm_borders(col = "red", lwd = 2) +
tm_view(set_zoom_limits = c(12,14))
<<<<<<< HEAD
=======
>>>>>>> 81336bfc1c4053bed01e1b7bc57e3309e34bd5e4